function [Y cache] = geodesic_best_new_cached(cache, state, reps)

% fast memory efficient implementation of geodesic interpolation algorithm which uses cached information
% state, the index of the element in the state sequence in that cache that is to be modified
% reps, the new number of repetition for the state

% restore some variables
clusters = cache.clusters;
SS = cache.SS;
repetitions = cache.repetitions;

% scalar variables
S = constructFullSequence(SS, repetitions);
T = sum(repetitions);
DIM = size(clusters.covars, 1) / 2;
Q = size(clusters.covars, 3);

% more variables
x = 1:DIM;
dx = DIM+x;

% restore more variables
A = cache.A;
B = cache.B;
C = cache.C;
D = cache.D;
mu = cache.mu;
dmu = cache.dmu;
AA = cache.AA;
BB = cache.BB;

% find out which part of the original cache can be reused.
reps_diff = reps - repetitions(state);
last_before_frame = sum(repetitions(1:(state-1))) + min([repetitions(state) reps]) + 1 - 1;
first_after_frame = sum(repetitions(1:state)) + 1 + 1;
range_to_compute = last_before_frame + (1:max([reps_diff+1 1]));

new_repetitions = repetitions;
new_repetitions(state) = reps;
newT = sum(new_repetitions);
newS = constructFullSequence(SS, new_repetitions);

start_index = 1;
recompute_start = (range_to_compute(1) == 1);
if recompute_start
    start_index = 2;
end

end_index = length(range_to_compute);
recompute_end = (range_to_compute(end) == (newT+1));
if recompute_end
    end_index = length(range_to_compute)-1;
end



% allocate memory for newAA
newAA = zeros(length(range_to_compute)*DIM, 3*DIM);
if recompute_start
    newAA(1:DIM, DIM+(1:2*DIM)) = [  D{newS(1)}, ...
                                   -(C{newS(1)} + D{newS(1)})];
end
for t = range_to_compute(start_index:end_index)-1
    rows = (t-range_to_compute(1)+1)*DIM + (1:DIM); cols = 1:3*DIM;
    newAA(rows, cols) = [-(B{newS(t)}   + D{newS(t)}), ...
                           A{newS(t)}   + B{newS(t)} + C{newS(t)} + D{newS(t)} + D{newS(t+1)}, ...
                         -(C{newS(t+1)} + D{newS(t+1)})];
end
if recompute_end
    rows = (newT-range_to_compute(1)+1)*DIM + (1:DIM);  cols = 1:2*DIM;
    newAA(rows, cols) = [-(B{newS(end)} + D{newS(end)}), ...
                           A{newS(end)} + B{newS(end)} + C{newS(end)} + D{newS(end)}];
end

before_rows = 1:(last_before_frame*DIM);
after_rows = ((first_after_frame-1)*DIM+1):size(AA, 1);
AAA = [AA(before_rows, :); newAA; AA(after_rows, :)];



% allocate memory for new BB
newBB = zeros(1, length(range_to_compute)*DIM);
if recompute_start
    newBB(1:DIM) = -C{newS(1)} * mu{newS(1)}(:) - D{newS(1)} * dmu{newS(1)}(:);
end
for t = range_to_compute(start_index:end_index)-1
    newBB((t-range_to_compute(1)+1)*DIM+(1:DIM)) =   (A{newS(t)}   + C{newS(t)}) * mu{newS(t)}(:) ...
                                                   + (B{newS(t)}   + D{newS(t)}) * dmu{newS(t)}(:) ...
                                                   -  C{newS(t+1)} * mu{newS(t+1)}(:) ...
                                                   -  D{newS(t+1)} * dmu{newS(t+1)}(:);
end
if recompute_end
    range = (newT-range_to_compute(1)+1)*DIM+(1:DIM);
    newBB(range) = (A{newS(end)} + C{newS(end)}) * mu{newS(end)}(:) + (B{newS(end)} + D{newS(end)}) * dmu{newS(end)}(:);
end
BBB = [BB(before_rows) newBB BB(after_rows)];



Y = reshape(solveBlockTridiagonal(AAA, BBB', DIM), DIM, newT+1)';



% setup cache
cache.repetitions = new_repetitions;
cache.AA = AAA;
cache.BB = BBB;
